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Theoretical  models  will  be  developed  for  thermomechanical  response 
in  terms  of  microstructural  parameters,  such  as  grain  size,  fiber 
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2.1  Overview 


Theoretical  work  during  the  past  year  was  concerned  with  (1) 
damage  growth  in  ceramic  matrices  at  high  temperature  and  with  (2) 
mixed-mode  delamination.  Improved  understanding  of  both  types  of 
damage  is  needed  in  order  to  predict  thermomechanical  behavior  of 
ceramic  composite  laminates. 


Under  item  (1),  models  of  crack  and  damage  growth  were 
developed  which  account  for  both  elastic  behavior  of  crystalline 
grains  and  viscous  behavior  of  grain  boundary  phases.  These  results 
for  viscoelastic  behavior  generalize  earlier  work  of  others  which 
appears  to  be  limited  to  the  extremes  of  elastic  behavior  and  viscous 
behavior.  This  analysis  has  not  yet  been  published,  and  therefore 
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essential  features  of  the  theory  and  some  results  are  given  in  the 
body  of  this  report  in  Section  2.2 

Item  (2)  on  delamination  was  written  and  published  during  the 
grant  period  [1];  the  publication  is  in  the  Appendix.  It  represents  a 
departure  from  other  formulations  of  mixed-mode  delamination 
analysis.  Traditionally,  the  total  energy  release  rate  and  its 
components  have  been  expressed  in  terms  of  moments  and  forces 
acting  on  the  boundaries  of  the  crack  (delamination)  tip  element  (cf. 
Fig.  3,  Appendix).  Instead,  we  developed  the  energy  release  rates 
directly  in  terms  of  crack  tip  moments  and  forces.  This  approach 
leads  to  simpler  results,  especially  for  complex  layups,  and  simplifies 
the  analytical  or  numerical  methods  needed  for  computing  mode 
ratios.  Some  further  work  on  the  prediction  of  mode  ratios  and  on 
the  effects  of  temperature  and  shear  deformations  is  currently 
underway. 

Experimental  studies  were  not  started  because  of  the  departure 
of  the  co-principal  investigator,  Professor  M.  C.  Lu,  from  Texas  A&M 
University  (shortly  after  the  grant  was  awarded)  and  Professor 
Schapery's  intention  to  move  to  The  University  of  Texas  at  Austin 
near  the  end  of  the  first  year  of  the  grant.  It  is  planned  to  initiate 
the  experimental  work  at  The  University  of  Texas  at  Austin  as  soon 
as  personnel  and  equipment  are  available. 
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2.2  Damage  Growth  in  Viscoelastic  Ceramics 

2.2.1.  Background.  An  important  micromechanism  for  time- 
dependent  deformation  in  many  ceramics  and  ceramic  composites  is 
the  nucleation,  growth  and  coalescence  of  cavities  and  crack-like 
defects  along  grain  boundaries  and  interfaces  [2].  Amorphous  phases 
that  are  often  present  along  grain  boundaries  and  interfaces  are 
usually  modeled  as  linear  or  nonlinear  viscous  materials  (i.e.,  the 
current  shear  stress  is  a  function  of  the  current  shear  strain  rate), 
and  we  shall  do  so  here.  The  objective  of  the  present  analysis  is  to 
develop  a  model  for  predicting  the  growth  of  individual  cracks  in 
viscous,  cavitating  grain  boundaries  or  interfaces  and  to  use  this 
model  to  predict  the  effect  of  a  distribution  of  such  cracks  on  the 
global  behavior  of  polycrystalline  ceramics. 

The  primary  background  references  for  this  work  are  the  crack 
tip  model  of  Thouless  et  al.  [3],  the  distributed  damage  model  of 
Suresh  and  Brockenbrough  [2],  and  the  viscoelastic  crack  growth 
theory  of  Schapery  [4],  In  [2]  the  effect  of  the  growth  of  distributed 
damage  on  creep  of  ceramics  and  ceramic  composites  is  predicted 
using  results  from  [3]  for  crack  growth  rate. 

Let  us  first  review  certain  aspects  of  these  two  important 
studies  before  discussing  the  new  theory.  Figure  1  shows  the  crack 
tip  model  of  [3],  in  which  the  cavities  grow  and  coalesce  in  a  very 
thin  viscous  layer  between  essentially  rigid  grains.  The  "viscous 
matrix"  outside  of  the  damage  zone  represents  the  undamaged, 
surrounding,  deformable  polycrystalline  material.  Grain 


deformations,  in  effect,  are  neglected  inside  and  outside  of  the 
damage  zone.  The  undamaged  continuum  deforms  through  straining 
of  the  viscous  grain  boundaries. 
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Fig.  1.-  The  theoretical  situation  considered 
in  [3].  The  damage  zone  consists  of  an  array 
of  grains  co-planar  with  the  crack  and 
embedded  in  a  linearly  viscous  matrix.  The 
growth  of  cavities  is  constrained  by  the 
viscous  matrix. 


The  result  for  crack  speed  &  may  be  written  obtained  in  the 

form 

a=^^G*(dA,ct/d,Af) 

where 


Ki 

— 

opening-mode  stress  intensity  factor 

d 

— 

grain  size 

a 

= 

length  of  damage  zone  (a  £  d) 

T1 

ss 

viscosity  of  grain  boundary  phase 

X 

= 

cavity  spacing  (d/X>l) 

Af 

critical  area  fraction  of  cavitated  grain 
boundary  at  which  cavities  coalesce  (for  the 
grain  at  the  crack  tip) 

G* 

— 

dimensionless  function 
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Graphs  of  the  theoretically  predicted  function  G*  are  given  in  [3]. 

It  is  of  particular  interest  to  observe  that  crack  speed  is 
proportional  to  K  and  VJ  and  inversely  proportional  to  r| .  This 
behavior  could  have  been  predicted  directly  from  dimensional 
analysis  (for  the  material  and  geometry  idealizations  assumed), 
although  this  approach  was  not  used  in  [3].  Rather,  a  detailed 
analysis  was  employed  which  not  only  provided  these  dependencies 
but  also  gave  the  function  G*  (graphically,  in  most  cases). 

Suresh  and  Brockenbrough  [2]  used  the  linear  relationship 
between  k  and  Ki  to  predict  the  creep  of  ceramics  when  all  of  the 
creep  is  due  to  distributed  microcrack  (damage)  growth.  They 
concluded  that  for  short  time  behavior  the  strain  rate  e  is 
proportional  to  Ki2  and  therefore 

e  ~  a2  (2) 

Prediction  and  experimental  results  for  polycrystalline  alumina  and 
alumina-silicon  carbide  composite  are  shown  in  Fig.  2;  their  theory 
gives  the  straight  lines  with  a  slope  of  two  in  the  figure. 

However,  in  [2]  the  ceramic  was  assumed  to  be  fully  elastic 
(except  for  fhe  viscous  layer  in  each  microcrack  tip  damage  zone) 
whereas  in  [3]  the  ceramic's  overall  behavior  was  assumed  to  be 
viscous  (with  rigid  grains).  Had  the  assumption  of  viscous  behavior 
been  used  in  [2],  the  lines  in  Fig.  2  would  have  had  a  slope  of  unity 
instead  of  two,  clearly  at  odds  with  the  experimental  data. 


Fig.  2.  Steady-state  creep  characteristics  of 
polycrystalline  alumina  and  alumina-silicon 
carbide  composite.  See  [2]  for  details.  The 
dark  dashed  lines  come  from  the  theory 
discussed  in  Section  2.2.4. 

Here,  we  shall  use  elastic  grains  and  a  viscous  grain  boundary 
phase  in  both  the  crack  tip  and  distributed  damage  models.  This 
generalization  will  enable  us  to  predict  behavior  which  is  in  better 
agreement  with  the  experimental  data. 

2.2.2.  The  Undamaged  Polvcrvstalline  Continuum.  Both  the 
crack  tip  model  and  the  model  for  global  material  response  with 
distributed  damage  require  a  constitutive  equation  for  the 
undamaged  portion  of  the  ceramic.  Linear  viscoelastic  behavior  will 
be  taken  into  account. 

For  a  uniaxial  stress,  the  axial  strain  for  an  arbitrary  stress 


history  o(x)  is. 


The  quantity  D(t)  is  the  "creep  compliance."  For  a  creep  test,  in 
which  the  stress  is  applied  at  i  =  0  and  held  constant  thereafter, 
equation  (3)  yields 

€  =  D(t)o  (4) 

Thus,  D(t)  may  be  obtained  from  creep  test  data  as  D  =  e/a  for  t  >  0; 

D  s  O  for  t  <  0.  Instead,  here  we  shall  predict  D(t)  as  well  as  the 
associated  creep  Poison's  ratio  vc(t),  using  the  micromechanical 
model  shown  in  Fig.  3. 

This  model  is  usually  called  the  "generalized  self-consistent 
scheme"  or  the  "three-phase  model"  in  the  context  of  particulate 
composites.  The  equations  for  the  overall  (or  effective)  elastic 
Young's  modulus  E  and  Poisson's  ratio  v  are  given  in  [5].  A  grain  is 
idealized  here  as  a  spherical  particle  and  the  grain  boundary  as  a 
concentric  matrix  layer.  A  representative  "two-phase  composite" 
consisting  of  a  grain  and  grain  boundary  phase  is  embedded  in  the 
third  phase;  the  latter  phase  has  the  properties  (E,  v)  the  undamaged 
ceramic.  These  properties  are  predicted  by  requiring  that  E  and  v  for 
the  three-phase  composite  in  Fig.  3  equal  those  for  an  effectively 
homogeneous  material  with  the  same  E  and  v.  The  equations  in  [5] 
may  be  used  to  calculate  these  elastic  constants  in  terms  of: 

Em*  Ep  =  matrix  and  particle  Young's  modulus, 
respectively. 

vm,  vp  =  matrix  and  particle  Poisson's  ratio, 
respectively. 

Vp  =  particle  volume  fraction 


Fig.  3.  Generalized  self-consistent  scheme 
without  microcracking.  The  particle  and 
matrix  represent  a  grain  and  grain  boundary 
phase,  respectively. 


In  turn,  the  correspondence  principle  of  linear  viscoelasticity  theory, 
along  with  Laplace  transform  inversion,  enables  us  to  predict  from 
these  elasticity  results  the  creep  functions  D(t)  and  vc(t)  for  a 

composite  with  viscoelastic  constituents.  Alternatively,  approximate 
results  for  D(t)  and  vc(t)  may  be  found  by  the  simpler  quasi-elastic 

method  [6],  which  we  use  here. 


We  assume  the  grains  are  elastic.  The  grain  boundary  phase 
(the  matrix)  is  assumed  elastic  in  dilatation  (constant  bulk  modulus) 
and  is  a  viscoelastic  Maxwell  model  in  shear, 

Jm  =  Jmo  +  tflm  (5) 

where 

Jm  =  shear  creep  compliance 

Jmo  =  initial  (elastic)  value  of  shear 

creep  compliance,  which  is  the 
reciprocal  of  the  initial  shear 
modulus. 

=  shear  viscosity 
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Observe  that  at  long  times  (t  »  JmoT|m)  the  grain  boundary  phase 
becomes  a  viscous  material,  while  at  short  times  it  is  elastic. 

Figures  4  and  5  are  graphs  of  the  logarithms  (base  10)  of  the 
plane-strain  creep  compliance  C(t)  of  the  polycrystalline  ceramic 
material  versus  time,  where 

T0  =  Jmotl,  the  creep  time  constant 
C  =  [1  -  vc2(t)]  D(t)  (6) 

and  C0  is  the  initial  value  of  C.  The  plane-strain  creep  compliance,  C 
(which  differs  only  slightly  from  the  uniaxial  creep  compliance,  D),  is 
shown  because  it  is  this  compliance  that  is  used  in  the  crack  growth 
theory  discussed  in  the  next  subsection. 

The  thickness  of  the  grain  boundary  phase  (b-a  in  Fig.  3) 
divided  by  the  grain  size  (the  mean  diameter  b+a  in  Fig.  3)  may  be 
expressed  in  terms  of  the  grain  (particle)  volume  fraction  vg  (=Vp). 
We  find  that  this  ratio,  denoted  by  h,  is 

h  =  (1  -  Vgl/3)/(l  +  vgl/3)  (7) 

For  the  values  of  vg  =  .9,  .99,  .999  (used  in  Figs.  4  and  5),  we  find, 
respectively,  h  =.018,  .0017,  .00017;  these  values  are  believed  to 
cover  a  realistic  range  for  ceramics  [2].  Similarly,  the  ratio  of  initial 
Young's  modulus  of  the  amorphous  grain  boundary  phase  to  that  for 
the  grains  in  actual  ceramics  is  likely  to  be  in  the  range  from  1/5  to 
unity;  these  are  the  values  used  for  Figs.  4  and  5  respectively. 


LOGCC/C0)  LOGCC/C0) 


Predictions  of  C/C0  as  a  function  of  initial  Poisson's  ratios  for 
both  phases  were  found  to  be  practically  independent  of  them  over 
the  range  from  zero  to  one-half;  the  value  of  0.3  was  used  for  both 
phases  in  these  figures.  As  noted  earlier,  the  grains  are  assumed 
elastic,  with  constant  values  of  Poisson's  ratio  and  Young's  modulus. 

At  very  long  times,  all  curves  in  Figs.  4  and  5  have  a  slope  of 
unity  (viz.  C~t)  as  expected,  since  the  deformation  is  dominated  by 
viscous  flow  of  the  grain  boundary  phase;  this  is  the  behavior 
assumed  in  the  crack  growth  theory  [3].  It  should  be  added  that  the 
long-time  shear  viscosity  tic,  say,  for  the  ceramic  may  be  found  using 
the  theory  in  [7], 

Tic  =  W32h3  (8) 


when  h«l. 

Especially  interesting  is  the  short  and  intermediate  time 
behavior  of  the  creep  compliance.  For  vg  =  .99  and  .999  there  is  a 
distinct  elastic  plateau  for  which  C/C0  =  2.3.  In  all  cases,  at  short 
times  there  is  a  nearly  constant  slope  region  in  which  slope  =  0.2, 
which  implies  C  ~  t  .  The  main  difference  between  the  results  in 
Figs.  4  and  5  is  a  time-shift.  Events  in  Fig.  4  occur  about  four  times 
faster  than  those  in  Fig.  5,  assuming  the  time  constant  T0  is  the  same 
for  all  cases. 

The  viscoelastic  behavior,  C  ~  t0-2,  and  elastic  behavior  where 
C/C0  ~  2.3,  stem  from  elastic  distortions  of  the  grains  as  they  interact 
with  the  grain  boundary  phase.  One  has  to  wait  a  relatively  very 
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long  time  until  C  ~  t,  when  the  grains  behave  as  relatively  rigid 
particles  in  a  viscous  matrix.  The  insensitivity  of  the  ceramic's  creep 
compliance  to  Poisson's  ratios  of  the  phases  implies  insensitivity  to 
the  bulk  moduli  of  the  grains  and  the  amorphous  grain  boundary 
phase. 

In  real  ceramics  the  thickness  of  the  grain  boundary  phase  is 
expected  to  vary  from  grain-to-grain.  This  implies  one  may  have  to 
account  for  a  distribution  of  local  volume  fractions  vg,  which  would 
lead  to  compliance  curves  that  are  combinations  of  those  in  Figs.  4  or 
5.  This  generalization  will  be  studied  in  future  work. 

2.2.3  Crack  Growth  Theory.  Equation  (1)  gives  the  crack  speed 
only  when  it  is  slow  enough  that  events  in  the  neighborhood  of  the 
crack  tip  take  place  in  the  very  long  time  regime,  C  ~  t.  We  may  now 
generalize  the  theory  to  allow  for  behavior  over  the  entire  time 
range  of  response  shown  in  Figs.  4  and  5  by  drawing  upon  Schapery's 
theory  in  [4]  using  the  idealized  crack  tip  model  in  Fig.  6. 


Figure  6.  Cross-section  of  an  idealized  crack 
used  in  [4],  The  failure  zone  here  is  to  be 
identified  with  the  damage  zone  in  Fig.  1  (i.c. 
the  cavities  in  a  viscous  layer  plus  the 
adjacent  rigid  grains. 


The  "failure  zone"  in  this  figure  is  to  be  identified  with  the  damage 
zone  of  Fig.  1,  which  consists  of  rigid  grains  and  a  cavitating  viscous 
layer.  The  theory  in  [4]  permits  us  to  use  mechanical  properties  for 
the  failure  zone  which  are  different  from  the  surrounding  ceramic 
continuum.  The  length  of  the  failure  zone  a  is  related  to  the  stress 
intensity  factor  Ki  according  to  [4,  Part  I], 

a  ~  K^/cfS  (9) 

where  Of  is  a  measure  of  the  normal  stress  acting  across  the  failure 
(or  damage)  zone.  From  [4,  Part  II]  the  crack  opening  displacement 
at  the  left  end  of  the  failure  zone  in  Fig.  6  is 

A  v  ~  C(ta)  K^/af  (10) 

where  C  is  the  plane-strain  creep  compliance  for  the  ceramic 
continuum  discussed  earlieT,  and 

t«  =  a/3d  (11) 

is  one-third  the  time  elapsed  when  the  crack  grows  an  amount  a. 

The  viscous  crack  growth  theory  in  [3]  uses  a  critical  cavity 
area  fraction  as  a  crack  growth  criterion;  according  to  [3,  Eq.  (5)],  this 
is  equivalent  to  a  critical  value  of  Av/X  and  thus  a  critical  value  Avc. 
From  Eqs.  (9)  -  (11), 

Avc~  C(a/3S)KiVoc  (12) 

For  a  viscous  body,  C  ~  a/3d,  and  if  a  is  constant  then  Eq.  (12) 
yields 

a-K1  (13) 

which  is  the  same  behavior  as  in  Eq.  (1). 
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For  a  power  law  viscoelastic  ceramic, 

C  =  (t/to)nC  i  (14) 

corresponding,  say,  to  the  short  time  region  in  Figs.  4  or  5  (for  which 
n  =  0.2)  and  if  a  is  constant,  then  Eq.  (12)  gives 

4~Ki1/n  (15) 


With  n  =  0.2, 

4  ~  Ki5  (16) 


It  was  assumed  in  [3]  and  in  developing  Eqs.  (13)-(16)  that  a  is 
constant.  However,  this  is  not  necessarily  correct.  Equation  (9)  must 
be  satisfied  (in  order  that  the  crack-tip  stress  not  be  infinite).  If  we 
suppose  that  all  geometric  parameters  at  the  crack  tip  are  constant 
(except  possibly  for  a)  then  for  a  viscous  failure  zone  af  ~  4,  which 

yields 

a  ~  Kj2/42  (17) 


Substitute  this  behavior  into  Eq.  (10)  after  using  Eq.  (14)  to  find 

4~KiP  (18) 


where 


P  = 


2n  +  2 
3n  +  1 


(19) 


For  n  =  1,  then  p  =  1  as  before;  this  result  and  Eq.  (17)  imply  a  is 
indeed  constant.  However,  if  n  *  1,  a  will  vary  with  Kr  For  an  elastic 


material,  n  =  0, 

4  ~  V  (20) 

whereas  if  n  =  0.2, 


4  -  Kji-5 


(21) 


which  gives  an  exponent  that  is  at  the  midpoint  between  elastic  and 
viscous  behavior!  It  should  be  noted  that  the  grain  boundary  phase 
is  viscous  in  [3]  but  viscoelastic  in  Eq.  (5)  and  Figs.  4  and  5.  For 
T/T0  >  10,  the  elastic  strain  is  negligible  and  therefore  the  two 
theories  are  for  viscous  behavior,  Jm  ~  t,  over  most  of  the  time  scale 
in  these  figures. 

A  more  detailed  analysis  leads  to  explicit  dependence  of  &  on 
the  various  physical  and  geometric  parameters  appearing  in  the 
cavitation  model  in  [3].  This  work  will  be  published  in  the  near 
future. 

2.2.4  Distributed  Damage  Model.  The  value  of  the  exponent  p 
has  a  significant  effect  on  the  overall  creep  strain  rate  of  a  ceramic 
with  growing  microcracks,  each  of  which  obeys  Eq.  (18).  In  order  to 
use  the  results  in  [2]  for  ceramics  which  are  elastic  (except  for  the 
effect  of  damage  growth)  we  suppose  the  linear  viscoelastic 
compliance  is  in  the  elastic  plateau  range  of  Figs.  4  and  5  (for  vg  > 
0.99,  say).  Then,  using  the  short-time  strain  rate  (due  to  damage) 
from  [2], 

e  ~  gp+1  (22) 

and  with  p  ~  1.5, 

€  ~  G2,5  (23) 

Whereas  if  even  shorter  effective  times  apply  to  the  crack  tip  region, 
we  use  elastic  behavior,  and  therefore  p  =  2  and 

e  ~  o3  (24) 

The  time  scale  for  the  compliance  around  the  damage  zone  is 


ta  ~  a/34,  which  is  smaller  than  the  time  scale  for  global  behavior.  It 
is  thus  appropriate  to  select  for  p  the  slope  at  times  less  than  those 
in  the  plateau. 

The  dark  dashed  lines  in  Fig.  2  correspond  to  Eq.  (22)  for  either 
elastic  behavior,  p  +  1  =  3,  or  viscoelastic  behavior,  p  +  1  =  2.5.  The 
experimental  data  do  not  all  have  the  same  slope,  which  may  be  due 
to,  at  least  in  part,  the  complexity  of  the  creep  compliance  of  a 
polycrystalline  ceramic  and  locally  varying  grain  volume  fraction  v  . 

o 

2.2.5  Conclusions.  These  results  on  crack  and  damage  growth 
are  very  encouraging  as  the  viscoelasticity  theory  leads  to  internally 
consistent  models  and  realistic  predictions,  all  expressed  in  terms  of 
basic  micromechanisms.  It  is  not  limited  to  constant  stress 
conditions.  Indeed,  the  theory  predicts  strong  frequency  effects  for 
crack  and  damage  growth  in  cyclic  loading,  such  as  reported  by  Han 
and  Suresh  [8].  Future  work  will  fill  in  details  of  this  approach  in 
order  to  relate  response  explicitly  to  micromechanical  parameters 
and  to  different  types  of  loading  and  thermal  histories.  Also,  the 
applicability  of  the  author's  earlier  work  on  distributed  damage 
growth  in  particulate  composites  [9]  and  a  low  temperature  ceramic 
(ice)  [10]  will  be  studied. 
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Prediction  of  the  energy  release  rate  (ERR)  and  Its  components  for  mixed¬ 
mode  delamination  of  composite  laminates  is  discussed.  A  classical  plate 
theory  (CPT)  version  of  Irwin's  virtual  crack  closure  method  is  developed 
and  used  for  the  ERR,  first  for  plane  strain  and  then  for  three- 
dimensional  deformations.  It  is  shown  that  CPT  does  not  provide  quite 
enough  information  to  obtain  a  decomposition  of  ERR  into  its  opening  and 
shearing  mode  components.  Results  from  a  continuum  analysis  are  needed  to 
complete  the  decomposition;  but  analysis  of  only  one  loading  case  is 
required  for  two-dimensional  and  certain  three-dimensional  problems.  In 
two  example  problems  the  finite  element  method  is  used  with  CPT  to 
complete  the  mode  decomposition.  Results  from  CPT  and  the  finite  element 
method  are  then  compared  for  several  cases. 


1.  INTRODUCTION 

Delamination  is  a  comnonly  occurring  type  of 
crack  growth  in  many  composite  laminates  .  The 
strength  and  stiffness  of  a  composite  structure 
may  be  significantly  affected  by  delamination, 
whether  the  loading  is  in-plane  or  out-of¬ 
plane.  One  indication  of  the  importance  of  the 
delamination  problem,  especially  for  fiber- 
reinforced  plastics,  is  the  large  number  of 
publications  devoted  to  the  subject,  such  as 
found  in  a  recent  book  on  fracture  of  composites 
(Friedrich,  1989).  The  energy  release  rate  (ERR) 
is  the  loading  parameter  commonly  used  to 
characterize  the  initiation  and  continuation  of 
crack  growth.  Delamination  growth  in  many 
composites  depends  on  not  only  the  (total)  ERR 
but  also  on  how  much  energy  is  in  the  shearing- 
and  opening-mode  components.  Numerous 

experimental  studies  (e.g.  Johnson  and 

Mangalglri,  1987)  have  shown  that  the  ERR 

required  for  delamination  growth  varys  with  the 

ratio  of  shearing-to-opening  mode  components  of 
the  ERR;  the  sensitivity  to  this  mode-ratio 
appears  to  be  least  in  the  toughest  composites. 
To  predict  initiation  of  growth,  one  compares  the 
available  and  required  values  of  ERR,  while 
stable  growth  rate  is  usually  expressed  in  terms 
of  the  available  ERR;  the  amplitude  of  ERR  Is 
often  used  for  cyclic  loading. 

In  this  paper  we  are  concerned  with  the 
prediction  of  ERR  and  Its  components  in  terms  of 
the  loading  acting  on  a  laminate.  For  use  with 
cyclic  loading  problems,  the  ERR  amplitude  and 

its  components  may  be  found  from  the  theory. 
There  is  no  fundamental  restriction  on  the 
geometry  of  the  laminate,  which  may  be  a  plate, 
shell  or  beam.  The  problem  may,  for  example,  be 
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that  of  predicting  edge  delamination  in  a 
uniaxial  tensile  specimen.  Fig.  1,  or  growth  of 
an  embedded  delamination  in  a  plate  under  in¬ 
plane  compression  or  a  transverse  load.  Fig.  2. 
Problems  like  the  latter  one  often  involve 
geometric  nonlinearities  due  to  the  out-of -plane 
deformations;  what  is  developed  here  can  be  used 
even  with  geometric  nonlinearities. 

The  approach  consists  of  analyzing  only  the 
neighborhood  of  a  short  segment  of  the 
delamination  edge;  the  cross-section  of  the  near¬ 
edge  region  is  illustrated  in  Fig.  3,  in  which 
there  is  one  delamination  or  crack  of  length  "a" 
within  this  element.  In  general,  the  delamina¬ 
tion  edge  is  a  curve  which  may  be  open  and  thus 
end  at  the  edges  of  the  laminate  (e.g.  Fig.  1)  or 
closed  (e.g.  Fig.  2).  The  x-axis  in  Fig.  3  is 
located  at  the  midsurface  and  is  normal  to  the 
local  delamination  edge.  We  shall  call  the 
geometry  in  Fig.  3  a  crack  tip  element.  The 
length  of  this  element  perpendicular  to  the  plane 
of  the  page  (the  y-direction)  is  denoted  by  L  ; 
it  is  assumed  to  be  short  enough  that  there  is  no 
significant  variation  in  the  loading  and  the 
edge's  orientation  with  respect  to  y.  It  is 
further  assumed  that  (a,b)  »  t,  but  yet  a  and  b 
are  small  enough  that  geometric  nonlinearities 
are  negligible,  given  the  moments  and  forces 
acting  on  this  element.  In  order  to  predict 
loading  on  the  element  (from  an  overall  plate  or 
shell  analysis)  it  may  be  necessary  to  account 
for  geometric  nonlinearities  and  possibly  other 
complicating  factors  such  as  shear  deforma¬ 
tions.  However,  for  the  crack  tip  element 
Itself,  we  assume  classical  linear,  elastic  plate 
theory  (CPT)  may  be  used  to  predict  deformations 
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FIG.  1.  Eight  ply  laminate  with  edge  delami¬ 
nations.  (After  Raju  et  al.,  1988). 


compression 

loads 


FIG.  2.  Locally  postbuckled  laminate  with  em¬ 
bedded  delamination.  (After  Whitcomb  and 
Shivakumar,  1989). 


and  strain  energies.  In  CPT  the  deformations  are 
defined  entirely  by  midsurface  strains  and 
curvatures;  external  forces  which  are  normal  to 
the  midsurface  and  act  at  locations  away  from 
this  element  are  fully  accounted  for  through  the 
loads  parallel  to  the  midsurface  and  the 
moments.  The  ERR  for  this  element  is  the  work 
(divided  by  L  aa)  that  becomes  available  at  the 
crack  tip  wherr  the  crack  advances  a  distance 
«a  along  the  length  Ly. 

Essentially  the*  same  considerations  have 
been  used  by  Storakers  and  Andersson  (1988), 
Whitcomb  (1989),  Whitcomb  and  Shivakumar  (1989) 
and  Williams  (1988)  to  derive  the  ERR.  In 
addition,  Williams  (1988)  decomposed  the  ERR  into 
opening  and  shearing  components;  however,  it  was 
assumed  that  a  pure  opening  mode  exists  when  the 
only  loading  is  =  -M2,  which  we  show  later  is 
not  generally  true.  In  order  to  determine  the 
components  of  ERR,  Whitcomb  (1984,  1986)  used 
finite  element  solutions  with  the  virtual  crack 
closure  method. 

In  this  paper  we  express  the  ERR  directly  in 
terms  of  the  one  concentrated  moment  and  two 
components  of  shearing  force  which  act  at  the 
delamination  edge.  These  three  reactions  are 
easily  written  in  terms  the  various  loads  and 
moments  which  act  on  the  boundaries  of  the  crack 
tip  element.  By  using  crack-tip  reactions,  we 
obtain  immediately  the  CPT  version  of  ERR 
components.  Also,  the  connection  between  these 
components  and  the  true  ERR  components  based  on  a 
local  continuum  analysis  of  the  crack  tip  is 
readily  established  in  terms  of  the  smallest 
possible  number  of  Independent  inputs.  To  the 
authors'  knowledge,  none  of  the  other 
publications  which  use  CPT  (Including  that  of 
Davidson  (1988),  which  is  a  forerunner  to  the 
present  analysis)  express  the  ERR  in  terms  of  the 
crack  tip  forces  and  moment.  Instead,  in  these 
other  works  the  ERR  is  written  in  terms  of  the 
several  forces  and  moments  which  act  on  the 
boundary  of  the  crack  tip  element. 

In  Section  2  the  ERR  is  derived  for  a  state 
of  plane  strain.  The  theory  is  extended  to 
three-dimensional  deformations  in  Section  3. 
Section  4  gives  examples  using  CPT  and 
corresponding  finite  element  predictions.  The 
finite  element  model  employed  is  described  in  the 
appendix. 

2.  PLANE  STRAIN  CRACK  TIP  ANALYSIS 


FIG.  3.  Crack  tip  element  for  plate  theory 
showing  loading  and  dimensions. 


Plate  Equations 

Assuming  a  state  of  plane  strain  for  now, 
and  referring  to  Fig.  3,  we  consider  only  the  in¬ 
plane  loads  N,  Nj  and  N2  and  moments  M,  Mj  and 
M2;  these  quantities  are  resultants  per  unit  of 
length  in  the  y-directlon,  and  are  positive  when 
they  act  in  the  directions  shown.  From  overall 
equilibrium  of  the  geometry  In  Fig.  3  and 
observing  that  zj  =  -t2/2  and  z2  =  tj/2, 

N  =  Nj+N2,  M  =  M1+M2+N2t1/2  -  Njt 2/2  (1) 

The  energy  release  rate  is  expressed  most 
simply  in  terms  of  the  crack  tip  shear  force  N- 
and  moment  M_,  Fig.  4.  The  presence  of  a  crack 
tip  is  fully  accounted  for  by  concentrated 
reactions  Nc  and  Mc  because,  in  the  context  of 
CPT,  there  are  no  tractions  which  act  across  the 
surface  of  crack  prolongation  in  the  uncracked 
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FIG.  4.  Free  body  diagram  for  plates  above 
and  below  the  crack  plane. 


section,  and  the  deformation  depends  only  upon 
moments  and  axial  forces.  Considering  the  top 
free-body  in  Fig.  4, 

Nc  =  -  Ni  +  Nj  (2) 


Mc  =  M1  +  "cV2  -  Mi 


(3) 


The  quantities  N,  and  M,  are  internal  stress- 
resultants  in  the^uncracked  section.  They  may  be 
found  from  CPT  and  expressed  in  terms  of  the 
overall  force  N  and  moment  M  in  this  section. 
Fig.  3,  as  follows.  Using  Jone's  (1975) 

notation,  for  the  full  laminate, 

N  =  Ae°+  B<  (4a) 

M  =  Be°+  Die  (4b) 

where  e°  is  the  midsurface  (z=0)  strain  and  <  is 
the  curvature.  Also,  A,B  and  D  are  the 
extensional  stiffness,  coupling  stiffness,  and 
bending  stiffness,  respectively.  Similarly,  for 
the  laminate  or  material  "1"  above  the  crack 
plane,  but  in  the  uncracked  section, 

N j-  Aie°u+  Bjk  =  AlE0+  (Br  ^t2/2)K  (5a) 

Mj=  BltJu+  DjK  =  Bie°+  (Or  B1t2/2)<  (5b) 

where  we  have  used  the  fact  that  the  midsurface 
strain  for  material  "1"  is 


0  0  0  i,  /o 

elu"  E  +  Z1K  -  E  -  Kt2/2 
How,  write  the  inverse  of  Eq.  (4)  as 

t°*  A'N  +  B'M 


(6) 

(7a) 


k  *  B'N  +  D'M 


(7b) 


in  which  the  primed  quantities  are  compliances. 
Use  this  result  in  Eq.  (5)  to  obtain  the  desired 
force  and  moment, 


Nr  anN  +  ai2M 

(8a) 

iN  +  *22^ 

(8b) 

where 

an=  AjA'  +  (Br  A1t2/2)B' 
a12=  AjB'  +  (Br  A1t2/2)D' 
a21=  BjA'  +  (Dr  81t2/2)B> 
a22=  B^'  +  (Dr  8^/2)  O' 

Equations  (2)  and  (3)  become 
Nc=  -H1+  auN  +  a12M 

Mc=  Mr  Njtj/2  +  (ajjtj/2  -  a21)H 

+  (a12tL/2  -  a22)M  (11) 


(9) 

(10) 


Energy  Release  Rate 

A  plate  theory  version  of  Irwin's  (1957) 
virtual  crack  closure  method  will  be  used  to 
derive  the  energy  release  rate.  First,  we 
consider  the  crack  tip  to  be  just  inside  the  left 
edge  of  the  geometry  in  Fig.  3.  The  left  edge  is 
then  the  right  boundary  of  the  uncracked  plate, 
and  may  be  treated  as  a  fixed  edge.  Now, 
referring  to  Fig.  4  with  Nc  and  Mc  initially 
zero,  gradually  increase  them  until  they  have  the 
values  needed  to  close  the  crack  along  the  length 
b;  these  values  are  given  by  Eqs.  (10)  and 
(11).  The  work  of  crack  closing  divided  by  b  is 
also  the  ERR,  G,  when  the  crack  length  grows  any 
small  amount  b.  Thus, 

G  =  2F  <Ncflu  +  Mcfie>  (12) 

where  au  =  u,-  u2  and  ae  =  e~-  e,.  The  u^  and 
0.  refer  to  Ttorfzontal  displacement  and  rotation 
of  the  crack  surfaces  for  the  top  (i=l)  and 
bottom  (i=2)  plates  due  only  to  Nc  and  Mc. 
Namely,  with  the  left  edge  fixed. 

Up  b(£j+K jt j/2) ,  9j=  Kjb  (13a) 

u2=  b(e2-K2t2/2) ,  e2=  <2b  (13b) 

where  e.  and  x.  are  midsurface  strain  and 
curvature1,  respectively,  due  to  Nc  and  Mc. 

For  the  material  above  the  crack  plane  the 
stiffness  matrix  has  components  Aj,  Bi  and  Dj.as 
in  Eq.  (5).  We  may  derive  the  midsurface  strain 
and  curvature  using  the  inverse  of  the  stiffness 
matrix  (A',  BS,  DM  and,  as  Inputs,  the 
midsurface  load  and  moment  shown  in  Fig.  5;  thus 


El=  AlNc+  Bi(Nctl/2  '  Mc)  (l4a) 

Kr  BiV  Di(NcV2  ■  Mc)  (14b) 

Substitution  of  Eq.  (14)  and  analogous  results 
for  the  material  below  the  crack  surface  into  Eq. 
(13)  yields,  finally, 

au/b  =  CjN^  c12Mc  (15a) 
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) 

) 


Ncti/2-Mc 


Nc  t2/2+Mc 


b 


and  the  are  real  coefficients.  If  the  are 
known,  the  G[  and  Gji  components  of  ERR  may  be 
calculated  in  terms  of  external  loads  and  moments 
through  Nc  and  M-,  Eqs.  (10)  and  (11). 

Thus,  it  only  remains  to  find  the  four  k^. 
By  equating  G  in  Eq.  (19)  to  Eg.  (17)  and 
matching  the  coefficients  of  N  ,  M  *  and  N_Mr  we 
find  c  c 


k2=  [c12kL+  s2(c1-k2)1/2c)/c1 

(21a) 

k3=  s3(cr  k2)1/2 

(21b) 

k4=  (c 32~  k^k2)/k3 

(21c) 

where  s2=  +1,  s3=  +1  and 

c  =  (c^  -  c22)1/2 

(21d) 

FIG.  5.  Loading  on  plates  above  and  below  the 
crack  plane  which  is  statically  equivalent  to 
loading  at  the  crack  tip. 


ae/b  =  c12Nc+  c2Mc  (15b) 

where 

cl=  AJ+  A£+  B^tr  B^t2+  Djt2/4  +  D£t2/4  (16a) 

c2=  D|+  D2  (16b) 

cl2=  D‘t2/2  -  D[t1/2  -  B‘-  B2  (16c) 

and  the  subscripts  1  and  2  on  (A',B',0')  refer  to 
the  plate  properties  above  and  below  the  crack 
plane,  respectively.  From  Eqs.  (12)  and  (15), 
the  ERR  becomes 

G  =  |  (clN2+  c/c  +  2c12NcMc)  (17) 

Given  that  G  >  0  for  Nc  and  Mc  not  both-zero,  it 
follows  that  c^  >  0,  c2  >  0  and  c^c2  >  c^. 


The  sign  S3  determines  the  sign  of  both  k3  and 
k*,  and  since  these  coefficients  define  the  sign 
of  Kjt  (and  g2)  we  may  use  S3  =  1  without  any 
real  Toss  in  generality.  For  the  cases  studied 
numerically  to-date 

c12^1  <  (cl-  ^  c  (22) 

and  thus  we  must  use  s2  =  1  so  that  k2  >  0; 
according  to  Eq.  (20a),  if  k2  <  0  a  sufficiently 
large  positive  crack  tip  moment  would  produce  the 
physically  unacceptable  result  of  a  negative 
opening-mode  stress  intensity  factor. 

There  is  not  enough  information  available 
from  CPT  to  find  the  one  remaining  unknown 
coefficient  kj.  As  kj  is  Independent  of  the 
loading,  we  may  use  numerical  or  other  results 
from  an  analysis  of  any  one  special  loading 
case.  It  is  helpful  to  select  applied  moments 
and  loads  which  make  N=M=0.  In  this  case  Eqs. 
(1),  (10)  and  (11)  yield 

N2=  -N1 ,  M2=  -M:+  Njt/2  (23) 

Nc=  -Np  Mc=  Mr  N1t1/2  (24) 


Mode  Decomposition 

For  those  cases  in  which  the  stresses  in  the 
continuum  in  the  neighborhood  of  the  crack  tip 
exhibit  the  behavior  r_i/  ,  where  r  is  distance 
from  the  tip,  we  may  write 


G  ■  Gj+  G j j 


(18) 


where  Gi  and  Gjt  are  the  mode  I  (opening  mode) 
and  mode  II  (shearing  mode)  energy  release 
rates.  In  turn  Gr~  K,  and  Grr-  K,f,  where  Kj  and 


IL 


are 


it-  r\ »■  ui urr  ,'[f i 

the  ‘essAciated11  stress  intensity 
factors.  These  factors  are  real,  linear 
functions  of  the  loads  which  produce  crack-tip 
stresses;  the  crack  tip  loading  here  is  defined 
entirely  by  Nc  and  Mc.  Thus,  we  may  write  Eq. 
(17)  in  the  form 


2G 


2  2 
*1+  92 


(19) 


gis  (2G j )**  =  kjNc+  k2Mc  (20a) 
g2=  (2Gj j ) **  -  k3Nc+  k4Mc  (20b) 


Further,  if  Ni  =  0  then  M2  =  -M,,  Nc  =  0  and 
M_  =  Mi.  Alternatively,  if  M,  =  N,ti/2  then  M  = 
0  and  M2  =  Njt2/2.  In  the  latter  case  Eq.  (20a) 
reduces  to  gi  =  -kjNj,  so  that  if  is  known  for 
this  one  loading  case, 

kp  -gj/Nj  (25) 
Nondimensional izatlon 

Here  we  shall  rewrite  the  primary  variables 
and  equations  in  dimensionless  form  to  allow  for 
thickness  and  modulus  scaling. 

Let  E  and  T  be  the  reference  modulus  and 
thickness,  respectively.  Then,  with  an  asterisk 
used  to  denote  a  nondimensionalized  variable,  we 
write  for  stiffnesses, 

A*  =  A/ET,  B*  =  B/ET2,  0*  =  D/ET3  (26) 
compliances, 

A‘ *  =  A'ET,  B'*  =  B‘ET2,  O'*  =  O'ET3 
load  and  moment. 


where 


(27) 
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N*  =  N/ET,  M*  =  M/ET2  (28) 

and  thickness  and  curvature, 

t*  =  t/T,  <*  =  kT  (29) 

where  the  subscript  "1"  or  "2"  may  be  used  with 
the  variables  In  Eqs.  (26)-(28),  as  appro¬ 
priate.  The  energy  release  rate  and  related 
variables  are 


G*  =  G/ET,  (30) 


-  CjET, 

c2*  =  c2ET3 

r  ♦  - 

,  c12 

c12£T' 

(31) 

9i*  = 

g./(ET)1/2 

for  i 

=  1,2 

(32) 

ki*  = 

k.j(ET)1/2 

for  i  = 

1,3 

(33) 

k.*  = 
Ki 

k,(ET)1/2T 

for  i  = 

2,4 

(34) 

The  ERR 

in  terms 

of  the 

dimensionless 

variables  is  the  same  as  that  in  Eq.  (17), 

G*  =  \  [c1*(Nc*)2+  c2*(Mc*)2+  2cx|  N*  Mc*1  (35) 

and  the  equations  for  the  coefficients  are  the 
same  as  in  Eq.  (16),  e.g. 

c2*  =  Dj*  +  D£*  (36) 

Similarly,  just  as  in  Eq.  (20), 

9l*  -  4*  Nc*  +  k2*  Mc*  (37) 

We  may  use  these  results  by  finding  solutions  for 
k*  using  E=T=1  and  then  generalizing  the  results 
to1  arbitrary  E  and  T  by  means  of  Eqs.  (33)  and 
(34);  e.g. 

kj-  k1*/(ET)1/2  (38) 

Homogeneous,  Isotropic  Plates 

In  the  numerical  examples  in  Section  4  the 
split  plate  is  assumed  to  be  made  of  one  or  two 
homogeneous,  isotropic  materials.  For  this  case 
the  stiffnesses  are  (Jones,  1975), 

At=  E1t1/(l-v2),  Bi=  0,  D.=  E-t?/12( 1-v? ) (39) 

where  v .  and  E.  are  Poisson's  ratio  and  Young's 
modulus,  respectively.  The  compliances  are 


a-  =  i/a1#  b;=  o,  o;  =  i/oi 

and  the  ERR  coefficients  are 

i  2  i  2 
1- V,  l-v« 

Cl-  t  +  £  t  I 
1  Vl  t2l2 
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are  obtained  by  omitting  the  terms. 

3.  THREE-DIMENSIONAL  CRACK  TIP  ANALYSIS 

We  continue  to  use  the  crack  tip  element  in 
Fig.  3,  but  now  allow  for  three  in-Dlanp  shear 
forces  In  the  y-direction  (say  NT,  NT,  Ny)  which 
act  on  the  same  three  areas  as  Nj,  and  N.  It 

is  assumed  that  Ly  »  b  »  t,  where  Ly  is  the 
element  length  in  the  y-direction,  but  that  L  is 
small  enough  for  all  loads  and  moments  ter  be 
essentially  independent  of  y  over  Ly. 

A  crack  tip  shear  force  in  the  y-direction, 
say  N'T,  also  may  exist.  Similar  to  Eq.  (2),  we 
obtain 

Nc  =  _N1  +  “l  (44) 

-  - 

In  order  to  find  NT  as  well  as  N.  and  M,  for  use 
in  Eqs.  (2)  and  (J),  the  full  three-dimensional 
plate  equations  are  needed.  In  place  of  Eqs.  (4) 
and  (5)  use  equations  of  the  type  (Jones,  1975), 


{51  -  [S !]{:'} 

(45) 

where  each  of  the  terms  N,M,e°,ic  represents  a  set 
of  three  components  (e.g.  Nx,  N  ,  N  ),  and  each 
element  A,B  and  D  represents  °a  3x3  symmetric 
matrix.  As  before,  one  set  of  equations  is 
needed  for  the  full  laminate  and  one  set  is 
needed  for  the  material  above  the  crack  plane. 

The  ERR  is  now 

G  =  2E  (Vu  +  V0  +  NcAV> 

(46) 

where  av  =  v,-  v2  is  the  difference 
displacements  across  the  crack  plane 
Mc  and  N^.  In  place  of  Eq.  (15),  use 

between  y- 
due  to  Nc, 

Au/b  =  CjNc+  c12Mc+  c13N> 

(47a) 

A9/b  =  C12Nc+  C2Mc+  C^N* 

(47b) 

Av/b  =  c13Nc+  c23Mc+  c3N^ 

(47c) 

where  the  c's  are  found  using  a  procedure  which 
is  the  same  as  for  the  plane  strain  case.  It 
should  be  noted  that  the  constraint  of  the 
uncracked  material  along  the  left  edge  of  the 
plates  in  Fig.  4  limits  the  deformation  caused  by 
Nc,  Mc  and  Ni  to  e°,  k  and  y°  ;  i.e.  the  c's 
are  to  be  round  *using  the  ’special  case 

V  V  Kxv=  0  in  Eq-  <45>- 
y  fronrtqs.  (46)  and  (47), 


(42) 

(43) 


Observe  that  c,?  *  0  when  the  top  and  bottom 
plates  are  identical.  Also,  plane  stress  results 


G  ■  i  Ic,n|-  c/s  c3!»»)?.  2c1;nc»c 


*  2c13NcN>,  2C„MC»»| 


When  there  Is  no  coupling  between  the 
force  N^  and  the  cross-sectional  variables 
au  and  a 8  (i.e.  Cjj  =  C23  =  0)  then  G  is  the  same 
as  in  Eq.  (17)  except  for  the  additional  term 
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c,(N^)  .  In  this  case,  and  without  (x.y) 
coupling  of  the  stress  singularities,  the  opening 
and  shearing  mode  components  (in  the  x-z  plane) 
of  ERR  depend  only  on  Nc  and  M-,  and  may  be  found 
from  Eq.  (20)  using  tne  coefficients  for  plane 
strain.  The  antipJane  shearing  component  of  G  is 
then  simply  c,(N“)V2. 

With  fulr  coupling  there  are  nine 


coefficients  in  the  three-dimensional  version 
of  Eq.  (20),  in  which  the  g^  (1=1 ,2 ,3)  are 
expressed  in  terms  of  N-,  Mj.,  and  N'.  Equation 
(48)  may  be  used,  as  in  the  plane  strain  case,  to 
solve  for  six  in  terms  of  the  remaining 
three.  However,  it  might  be  easier  to  use  a 
method  similar  to  Whitcomb's  (1984)  approach; 
three  different  loading  cases  (to  obtain  three 
different  sets  of  Nc,  Mc  and  N-)  would  be 
analyzed  by  the  finite  element  method  to  find  the 
nine  k^. 


4.  EXAMPLES 


Two-Layer  Plate 

Consider  a  crack  between  two  different 
homogeneous  isotropic  plates  of  equal  thickness, 
for  which  we  select  tj  =  t2  =  1,  E^  =  1,  E2  = 
1.923,  \>,=  0.3,  v-=  0.  For  plane  strain  and 

these  values  of  Poisson's  ratio  and  modulus  the 
stress  singularity  at  the  crack  tip  is  real  (Raju 
et  al.,  1988)  and  thus  does  not  have  the 
oscillatory  character  common  to  cracks  between 
dissimilar  media. 

Equations  (41 )-(43)  provide  the  coefficients 
in  the  ERR,  Eq.  (17);  we  find  that 

G  =  |(5.72  N^+  17.2  4.68  N^)  (49) 

The  coefficient  k,  in  Eq.  (20a)  is  found  to 
be  -0.174  from  finite  element  analysis  (cf. 
Appendix)  and  Eq.  (25).  The  loading  used  for 
this  case  is  Nj  =  2,  N2  =  -2,  Mj  =  M2  =  1,  which 
yields  Nc  =  -2  and  M  =  0.  The  remaining  ki  are 
obtained  from  Eq.  (21)  to  give  in  Eq.  (20), 

g1=  ( 2G j ) ^  -0.174Nc+  4.09Mc  (50a) 

g2=  ( 2G j j )*5=  2.39Nc-  0.683MC  (50b) 

Table  I  compares  the  plate  predictions  (P) 
from  Eqs.  (49)  and  (50)  with  values  determined 
directly  from  the  finite  element  model  (F)  for 
several  different  loading  cases  in  which  Nj,  N2, 
Mj  and  M2  are  varied. 


TABLE  I.  ERR  COMPARISONS  FOR  THE  TWO¬ 
LAYERED  PLATE 


Case 

Nc 

Mc 

Gp 

Gf 

6I 

GII 
u  p 

ci 

G"r 

1 

-2 

0 

11.4 

11.5 

5.34x10" 

2 

-3.09 

.152 

28.6 

28.6 

.0240 

.0240 

3 

-6.83 

1.80 

190 

191 

.237 

.241 

4 

-1.78 

.712 

16.3 

16.4 

.465 

.464 

5 

-4.39 

3.22 

177 

177 

1.21 

1.21 

6 

-3.81 

3.84 

203 

202 

1.95 

1.94 

7 

-1.06 

1.36 

22.3 

22.2 

2.75 

2.73 

8 

0 

1 

8.58 

8.57 

35.8 

36.0 

It  is  seen  that  agreement  between  plate  theory 
and  finite  element  analysis  is  very  good  for  all 
cases. 

Homogeneous  Plate 

The  second  problem  is  that  of  a  homogeneous, 
isotropic  plate  in  plane  strain  with  an  off- 
center  crack.  Specifically,  ti  =  1,  t2  =  0.5,  Ei 
=  Eo  =  1,  v,=  v,=  0.3.  Use  of  Eqs.  (4I)-(43)  ana 
(17f  yields,1  £ 

G  =  |(10.9N*+  98.3M*+  32.8NCMC)  (51) 

Equation  (25)  and  a  finite  element  analysis  give 
kj^  =  0.398,  so  that  from  Eqs.  (20)  and  (21), 

g1=  0.398NC+  9.12MC  (52a) 

g2=  3.28NC+  3.89MC  (52b) 

The  kj  value  was  found  from  the  pure  shear 
loading  case  Nc  =  2,  Mc  =  0.  One  important 
comparison  that  can  be  made  between  plate  theory 
and  finite  element  analysis  is  for  the  other 
extreme  case  of  pure  moment  loading.  With  N_  = 
0,  M-  =  1  we  found  very  good  agreement,  in  that 
both  methods  yielded  G  =  49.1  and  Gj/Gjj  = 
5.50.  Other  loading  cases  were  analyzed,  with 
the  same  conclusion. 

Discussion 

The  many  loading  cases  studied  numerically 
in  both  problems  have  used  a  variety  of  values 
for  N^,  N2,  M,  and  M2.  That  the  only  effective 
crack  tip  loading  parameters  are  N_  and  Mc  was 
confirmed  by  the  agreement  of  tne  continuum 
analysis  (using  finite  elements)  with  CPT.  That 
CPT  provides  total  ERR  which  is  usually  in  good 
agreement  with  finite  element  predictions  has 
also  been  reported  by  Whitcomb  (1989)  for  a 
postbuckled,  embedded  delamination.  The  accuracy 
of  CPT  in  predicting  ERR  thus  makes  it  a  useful 
tool  for  checking  whether  or  not  the  mesh  size 
and  proportions  in  finite  element  models  are 
acceptable. 

In  the  two  example  problems  the  values  tj  = 
Ej  =  1  were  used.  By  selecting  T  =  tj  and  E  =  Ej 
in  the  earlier  section  on  nondimensionalization 
we  may  convert  the  calculated  k^  and  c^  values, 
which  may  be  interpreted  as  kt  and  ct,  to  those 
for  other  moduli  and  thicknesses  using,  for 
example,  Eq.  (38). 

It  is  of  Interest  to  observe  that  Eqs.  (50a) 
and  (52a)  show  Gj  does  not  necessariliy  vanish 
even  if  the  crack  tip  moment  vanishes. 
Similarly,  when  the  shear  force  vanishes  G^  is 
not  necessarily  zero.  Such  results  are  possible 
because  N.  and  M.  are  only  resultants,  and  thus 
do  not  reflect  details  in  the  local  stress 
distribution.  On  the  other  hand  kj  is  relatively 
small  in  both  problems,  and  it  would  be 
interesting  to  investigate  whether  or  not  k^  =  0 
can  be  used  generally.  If  this  approximation  is 
valid,  then  k2,  k3  and  k4  may  be  derived  directly 
from  Eq.  (21)7  without  having  to  use  results  from 
finite-element  or  other  analyses.  The  accuracy 
of  delamination  predictions  with  kj  =  0  depends 
on  the  delamination  criterion.  For  example,  if 
the  critical  value  of  G  is  independent  of  mode 
ratio,  then  any  value  of  kj  <  c^  will  suffice 


» 
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because  Eq.  (21)  guarantees  that  6,  Eq.  (17),  is 
independent  of  kj.  Note  also  that  symmetry 
yields  kj  =  =  ci2  =  0  exactly  for  a  symmetric 

laminate  with  a  midsurface  crack. 

The  use  of  kj  =  0  avoids  the  problem  of 
finding  Gj/Gtj  from  a  continuum  analysis  when  the 
stress  singularity  is  oscillatory.  Raju  et  al. 
(1988)  have  shown  that  Gj  and  Gjt  do  not  converge 
(although  G  does  converge)  with  decreasing  finite 
element  mesh  size;  an  analytically  evaluated 
limit  shows  similar  behavior.  Thus,  kj  cannot  be 
found  for  this  common  type  of  singularity  at  a 
crack  tip  between  dissimilar  media.  Further 
study  is  needed  to  determine  if  the  use  of  k^  =  0 
or  any  other  real  value  is  a  valid  way  of 
avoiding  the  oscillatory  singularity  problem. 

5.  CONCLUSIONS 

The  energy  release  rate  (ERR)  has  been 
derived  for  two-  and  three-dimensional  delami¬ 
nation  problems  using  classical  plate  theory 
(CPT).  The  plate  analysis  was  simplified  by 
recognizing  that  the  three-dimensional  ERR  and 
its  components  can  be  fully  determined  from  one 
moment  and  two  shear  forces  at  the  crack  tip;  for 
two  dimensions,  there  is  only  one  moment  and  one 
shear  force.  Decomposition  of  the  plate  theory 
ERR  into  opening  and  shearing  mode  components  was 
accomplished  by  combining  it  with  numerical 
results  from  a  continuum  analysis;  an  approxima¬ 
tion  which  avoids  the  use  of  any  solutions  beyond 
plate  theory  was  proposed.  Finally,  good 
agreement  between  CPT  and  finite  element 
predictions  was  shown  for  ERR  and  its  components. 
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APPENDIX 


Figures  6  and  7  show  the  finite  element 
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FIG.  7.  Expanded  view  of  crack  tip  neighbor¬ 
hood  . 


model  employed  to  obtain  results  for  plane  strain 
given  in  Section  4.  We  used  four  node  isopara¬ 
metric,  quadrilateral  elements  in  the  code 
MSC»PAL2  Version  3.5  (The  MacNeal-Schwendler 

Corp.,  Los  Angeles,  CA  90041)  and  pre-  and  post¬ 
processor  FEMAP  Version  3.0  (Enterprise  Software 
Products,  Inc.,  Harleysville,  Pa.,  19438).  The 
parameters  tj  and  E^  were  selected  as  the 

reference  thickness  and  reference  modulus, 

respectively,  and  thus  we  used  for  calculation 
purposes  tj  =  E1  =  1.  Also,  a=b=10.  For  the 

case  t2  =  .5,  the  same  mesh  as  in  Figs.  6  and  7 
was  used  but  the  three  bottom  rows  of  elements 
were  removed.  Loads  and  moments  were  applied  to 
the  right  end  of  the  model  by  specifying  values 
for  four  loads  at  the  top  and  bottom  nodes  on  the 
plates  above  and  below  the  crack.  Crack  tip 
nodel  forces  were  found  by  using  very  stiff 
horizontal  and  vertical  spring  elements  to 
connect  nodes  across  the  crack  plane;  the 
displacement  between  these  nodes  was  at  least 
four  orders  of  magnitude  less  than  that  for  the 
adjacent  crack  surface  nodes. 

In  order  to  determine  if  the  mesh  was  fine 
enough  to  provide  accurate  results,  additional 
runs  were  made  using  elements  around  the  crack 
tip  of  size  0.01x0.01  in  place  of  the  0.05x0.05 
elements  shown  in  Fig.  7.  In  most  cases  the 
change  was  less  than  one  percent  for  G  and  five 
percent  for  Gj  and  Gjj. 


